metaSEM package
OpenMx 2.11.5, and metaSEM 1.2.0.## Required packages for this workshop
lib2install <- c("metaSEM", "semPlot", "readxl")
## Install them automatically if they are not available on your computer
for (i in lib2install) {
if (!(i %in% rownames(installed.packages()))) install.packages(i)
}
“for their development of models, statistical procedures, and a computer algorithm (LISREL by Karl G. Joreskog; EQS by Peter M. Bentler) for Structural Equation Modeling (SEM) that changed the way in which inferences are made from observational data by permitting hypotheses derived from theory to be tested.”
## Asymptotic sampling covariance matrix based on the harmonic mean (375):
## x2_x1 x3_x1 x3_x2
## x2_x1 0.0010923 0.0004864 0.0004864
## x3_x1 0.0004864 0.0010923 0.0004864
## x3_x2 0.0004864 0.0004864 0.0010923
## Asymptotic sampling covariance matrix calculated in TSSEM approach:
## x2_x1 x3_x1 x3_x2
## x2_x1 0.0020480 0.0005768 0.0004079
## x3_x1 0.0005768 0.0008192 0.0002580
## x3_x2 0.0004079 0.0002580 0.0004096
metaSEM package provides two methods to achieve this:
diag.constraints=FALSE);diag.constraints=TRUE)RE.type="Diag").metaSEM package.
tssem1()tssem2()tssem2() automatically handles whether a fixed- or random-effects model is used in the stage 1 analysis.metaSEM package also provides a function lavaan2RAM() to convert the lavaan syntax into RAM specification.metaSEM package as Digman97.Info) includes the study characteristics whereas the other sheets (Study 1 to Study 14) store the correlation matrices.## Load the library to read XLSX file
library(readxl)
## Read the study characteristics
study <- read_xlsx("Digman97.xlsx", sheet="Info")
head(study)
## # A tibble: 6 x 3
## Study n Cluster
## <chr> <dbl> <chr>
## 1 Digman 1 (1994) 102 Children
## 2 Digman 2 (1994) 149 Children
## 3 Digman 3 (1963c) 334 Children
## 4 Digman & Takemoto-Chock (1981b) 162 Children
## 5 Graziano & Ward (1992) 91 Adolescents
## 6 Yik & Bond (1993) 656 Young adults
## Create an empty list to store the correlation matrices
Digman97.data <- list()
## Read 1 to 14 correlation matrices
for (i in 1:14) {
## Read each sheet and convert it into a matrix
mat <- as.matrix(read_xlsx("Digman97.xlsx", sheet=paste0("Study ", i)))
## Add the row names
rownames(mat) <- colnames(mat)
## Save it into a list
Digman97.data[[i]] <- mat
}
## Add the names of the studies
names(Digman97.data) <- study$Study
## Show the first few studies
head(Digman97.data)
## $`Digman 1 (1994)`
## A C ES E I
## A 1.00 0.62 0.41 -0.48 0.00
## C 0.62 1.00 0.59 -0.10 0.35
## ES 0.41 0.59 1.00 0.27 0.41
## E -0.48 -0.10 0.27 1.00 0.37
## I 0.00 0.35 0.41 0.37 1.00
##
## $`Digman 2 (1994)`
## A C ES E I
## A 1.00 0.39 0.53 -0.30 -0.05
## C 0.39 1.00 0.59 0.07 0.44
## ES 0.53 0.59 1.00 0.09 0.22
## E -0.30 0.07 0.09 1.00 0.45
## I -0.05 0.44 0.22 0.45 1.00
##
## $`Digman 3 (1963c)`
## A C ES E I
## A 1.00 0.65 0.35 0.25 0.14
## C 0.65 1.00 0.37 -0.10 0.33
## ES 0.35 0.37 1.00 0.24 0.41
## E 0.25 -0.10 0.24 1.00 0.41
## I 0.14 0.33 0.41 0.41 1.00
##
## $`Digman & Takemoto-Chock (1981b)`
## A C ES E I
## A 1.00 0.65 0.70 -0.26 -0.03
## C 0.65 1.00 0.71 -0.16 0.24
## ES 0.70 0.71 1.00 0.01 0.11
## E -0.26 -0.16 0.01 1.00 0.66
## I -0.03 0.24 0.11 0.66 1.00
##
## $`Graziano & Ward (1992)`
## A C ES E I
## A 1.00 0.64 0.35 0.29 0.22
## C 0.64 1.00 0.27 0.16 0.22
## ES 0.35 0.27 1.00 0.32 0.36
## E 0.29 0.16 0.32 1.00 0.53
## I 0.22 0.22 0.36 0.53 1.00
##
## $`Yik & Bond (1993)`
## A C ES E I
## A 1.00 0.66 0.57 0.35 0.38
## C 0.66 1.00 0.45 0.20 0.31
## ES 0.57 0.45 1.00 0.49 0.31
## E 0.35 0.20 0.49 1.00 0.59
## I 0.38 0.31 0.31 0.59 1.00
## Extract the sample sizes
Digman97.n <- study$n
Digman97.n
## [1] 102 149 334 162 91 656 70 70 277 227 1000 227 91 1040
## Extract the cluster
Digman97.cluster <- study$Cluster
Digman97.cluster
## [1] "Children" "Children" "Children" "Children" "Adolescents" "Young adults" "Young adults" "Young adults" "Mature adults"
## [10] "Mature adults" "Mature adults" "Mature adults" "Mature adults" "Mature adults"
0 and 1 means that the values are fixed at 0 and 1, respectively."0.3*Alpha_A": the factor loading from Alpha to A is free with a starting value of 0.3.model <- "## Factor loadings
## Alpha is measured by A, C, and ES
Alpha =~ A + C + ES
## Beta is measured by E and I
Beta =~ E + I
## Factor correlation between Alpha and Beta
Alpha ~~ Beta"
## Display the model
plot(model, color="yellow")
## Convert the lavaan syntax into a RAM model as the metaSEM only knows the RAM model
## It is important to ensure that the variables are arranged in A, C, ES, E, and I.
RAM <- lavaan2RAM(model, obs.variables=c("A","C","ES","E","I"),
A.notation="on", S.notation="with")
RAM
## $A
## A C ES E I Alpha Beta
## A "0" "0" "0" "0" "0" "0*AonAlpha" "0"
## C "0" "0" "0" "0" "0" "0*ConAlpha" "0"
## ES "0" "0" "0" "0" "0" "0*ESonAlpha" "0"
## E "0" "0" "0" "0" "0" "0" "0*EonBeta"
## I "0" "0" "0" "0" "0" "0" "0*IonBeta"
## Alpha "0" "0" "0" "0" "0" "0" "0"
## Beta "0" "0" "0" "0" "0" "0" "0"
##
## $S
## A C ES E I Alpha Beta
## A "0*AwithA" "0" "0" "0" "0" "0" "0"
## C "0" "0*CwithC" "0" "0" "0" "0" "0"
## ES "0" "0" "0*ESwithES" "0" "0" "0" "0"
## E "0" "0" "0" "0*EwithE" "0" "0" "0"
## I "0" "0" "0" "0" "0*IwithI" "0" "0"
## Alpha "0" "0" "0" "0" "0" "1" "0*AlphawithBeta"
## Beta "0" "0" "0" "0" "0" "0*AlphawithBeta" "1"
##
## $F
## A C ES E I Alpha Beta
## A 1 0 0 0 0 0 0
## C 0 1 0 0 0 0 0
## ES 0 0 1 0 0 0 0
## E 0 0 0 1 0 0 0
## I 0 0 0 0 1 0 0
##
## $M
## A C ES E I Alpha Beta
## 1 0 0 0 0 0 0 0
## method="FEM": fixed-effects TSSEM
fixed1 <- tssem1(Digman97.data, Digman97.n, method="FEM")
## summary of the findings
summary(fixed1)
##
## Call:
## tssem1FEM(Cov = Cov, n = n, cor.analysis = cor.analysis, model.name = model.name,
## cluster = cluster, suppressWarnings = suppressWarnings, silent = silent,
## run = run)
##
## Coefficients:
## Estimate Std.Error z value Pr(>|z|)
## S[1,2] 0.363278 0.013368 27.1760 < 2.2e-16 ***
## S[1,3] 0.390375 0.012880 30.3076 < 2.2e-16 ***
## S[1,4] 0.103572 0.015048 6.8830 5.862e-12 ***
## S[1,5] 0.092286 0.015047 6.1331 8.621e-10 ***
## S[2,3] 0.416070 0.012519 33.2344 < 2.2e-16 ***
## S[2,4] 0.135148 0.014776 9.1464 < 2.2e-16 ***
## S[2,5] 0.141431 0.014866 9.5135 < 2.2e-16 ***
## S[3,4] 0.244459 0.014153 17.2724 < 2.2e-16 ***
## S[3,5] 0.138339 0.014834 9.3259 < 2.2e-16 ***
## S[4,5] 0.424566 0.012375 34.3070 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Goodness-of-fit indices:
## Value
## Sample size 4496.0000
## Chi-square of target model 1505.4443
## DF of target model 130.0000
## p value of target model 0.0000
## Chi-square of independence model 4471.4242
## DF of independence model 140.0000
## RMSEA 0.1815
## RMSEA lower 95% CI 0.1736
## RMSEA upper 95% CI 0.1901
## SRMR 0.1621
## TLI 0.6580
## CFI 0.6824
## AIC 1245.4443
## BIC 412.0217
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## extract coefficients
coef(fixed1)
## A C ES E I
## A 1.00000000 0.3632782 0.3903748 0.1035716 0.09228557
## C 0.36327824 1.0000000 0.4160695 0.1351477 0.14143058
## ES 0.39037483 0.4160695 1.0000000 0.2444593 0.13833895
## E 0.10357155 0.1351477 0.2444593 1.0000000 0.42456626
## I 0.09228557 0.1414306 0.1383390 0.4245663 1.00000000
fixed2 <- tssem2(fixed1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)
summary(fixed2)
##
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj),
## n = sum(tssem1.obj$n), Amatrix = Amatrix, Smatrix = Smatrix,
## Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis,
## intervals.type = intervals.type, mx.algebras = mx.algebras,
## model.name = model.name, suppressWarnings = suppressWarnings,
## silent = silent, run = run)
##
## 95% confidence intervals: z statistic approximation
## Coefficients:
## Estimate Std.Error lbound ubound z value Pr(>|z|)
## AonAlpha 0.562800 0.015380 0.532655 0.592944 36.593 < 2.2e-16 ***
## ConAlpha 0.605217 0.015324 0.575183 0.635251 39.495 < 2.2e-16 ***
## EonBeta 0.781413 0.034244 0.714297 0.848529 22.819 < 2.2e-16 ***
## ESonAlpha 0.719201 0.015685 0.688458 0.749943 45.852 < 2.2e-16 ***
## IonBeta 0.551374 0.026031 0.500354 0.602394 21.181 < 2.2e-16 ***
## AlphawithBeta 0.362678 0.022384 0.318806 0.406550 16.203 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Goodness-of-fit indices:
## Value
## Sample size 4496.0000
## Chi-square of target model 65.4526
## DF of target model 4.0000
## p value of target model 0.0000
## Number of constraints imposed on "Smatrix" 0.0000
## DF manually adjusted 0.0000
## Chi-square of independence model 3112.7380
## DF of independence model 10.0000
## RMSEA 0.0585
## RMSEA lower 95% CI 0.0465
## RMSEA upper 95% CI 0.0713
## SRMR 0.0284
## TLI 0.9505
## CFI 0.9802
## AIC 57.4526
## BIC 31.8088
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
semPlot package (Epskamp, 2014).24plot(fixed2, color="green")
# Display the original study characteristic
table(Digman97.cluster)
## Digman97.cluster
## Adolescents Children Mature adults Young adults
## 1 4 6 3
## Younger participants: "Children" and "Adolescents"
## Older participants: "Mature adults"
sample <- ifelse(Digman97.cluster %in% c("Children", "Adolescents"),
yes="Younger participants", no="Older participants")
table(sample)
## sample
## Older participants Younger participants
## 9 5
## cluster: variable for the analysis with cluster
fixed1.cluster <- tssem1(Digman97.data, Digman97.n, method="FEM", cluster=sample)
summary(fixed1.cluster)
## $`Older participants`
##
## Call:
## tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis,
## model.name = model.name, suppressWarnings = suppressWarnings)
##
## Coefficients:
## Estimate Std.Error z value Pr(>|z|)
## S[1,2] 0.297471 0.015436 19.2716 < 2.2e-16 ***
## S[1,3] 0.370248 0.014532 25.4774 < 2.2e-16 ***
## S[1,4] 0.137694 0.016403 8.3944 < 2.2e-16 ***
## S[1,5] 0.098061 0.016724 5.8637 4.528e-09 ***
## S[2,3] 0.393692 0.014146 27.8306 < 2.2e-16 ***
## S[2,4] 0.183045 0.016055 11.4009 < 2.2e-16 ***
## S[2,5] 0.092774 0.016643 5.5743 2.485e-08 ***
## S[3,4] 0.260753 0.015554 16.7645 < 2.2e-16 ***
## S[3,5] 0.096141 0.016573 5.8009 6.596e-09 ***
## S[4,5] 0.411756 0.013900 29.6225 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Goodness-of-fit indices:
## Value
## Sample size 3658.0000
## Chi-square of target model 825.9826
## DF of target model 80.0000
## p value of target model 0.0000
## Chi-square of independence model 3000.9731
## DF of independence model 90.0000
## RMSEA 0.1515
## RMSEA lower 95% CI 0.1424
## RMSEA upper 95% CI 0.1611
## SRMR 0.1459
## TLI 0.7117
## CFI 0.7437
## AIC 665.9826
## BIC 169.6088
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
##
## $`Younger participants`
##
## Call:
## tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis,
## model.name = model.name, suppressWarnings = suppressWarnings)
##
## Coefficients:
## Estimate Std.Error z value Pr(>|z|)
## S[1,2] 0.604330 0.022125 27.3141 < 2.2e-16 ***
## S[1,3] 0.465536 0.027493 16.9327 < 2.2e-16 ***
## S[1,4] -0.031381 0.035940 -0.8732 0.38258
## S[1,5] 0.061508 0.034547 1.7804 0.07500 .
## S[2,3] 0.501432 0.026348 19.0311 < 2.2e-16 ***
## S[2,4] -0.060678 0.034557 -1.7559 0.07911 .
## S[2,5] 0.320987 0.031064 10.3330 < 2.2e-16 ***
## S[3,4] 0.175437 0.033675 5.2097 1.891e-07 ***
## S[3,5] 0.305149 0.031586 9.6609 < 2.2e-16 ***
## S[4,5] 0.478640 0.026883 17.8045 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Goodness-of-fit indices:
## Value
## Sample size 838.0000
## Chi-square of target model 346.2810
## DF of target model 40.0000
## p value of target model 0.0000
## Chi-square of independence model 1470.4511
## DF of independence model 50.0000
## RMSEA 0.2139
## RMSEA lower 95% CI 0.1939
## RMSEA upper 95% CI 0.2355
## SRMR 0.1411
## TLI 0.7305
## CFI 0.7844
## AIC 266.2810
## BIC 77.0402
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
Beta_I=3.29, for the younger participants are problematic.fixed2.cluster <- tssem2(fixed1.cluster, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)
summary(fixed2.cluster)
## $`Older participants`
##
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj),
## n = sum(tssem1.obj$n), Amatrix = Amatrix, Smatrix = Smatrix,
## Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis,
## intervals.type = intervals.type, mx.algebras = mx.algebras,
## model.name = model.name, suppressWarnings = suppressWarnings,
## silent = silent, run = run)
##
## 95% confidence intervals: z statistic approximation
## Coefficients:
## Estimate Std.Error lbound ubound z value Pr(>|z|)
## AonAlpha 0.512656 0.018206 0.476973 0.548340 28.158 < 2.2e-16 ***
## ConAlpha 0.549967 0.017945 0.514795 0.585140 30.647 < 2.2e-16 ***
## EonBeta 0.967136 0.058751 0.851986 1.082287 16.462 < 2.2e-16 ***
## ESonAlpha 0.732174 0.018929 0.695074 0.769273 38.680 < 2.2e-16 ***
## IonBeta 0.430649 0.029634 0.372568 0.488730 14.532 < 2.2e-16 ***
## AlphawithBeta 0.349236 0.028118 0.294125 0.404346 12.420 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Goodness-of-fit indices:
## Value
## Sample size 3658.0000
## Chi-square of target model 21.9954
## DF of target model 4.0000
## p value of target model 0.0002
## Number of constraints imposed on "Smatrix" 0.0000
## DF manually adjusted 0.0000
## Chi-square of independence model 2273.3236
## DF of independence model 10.0000
## RMSEA 0.0351
## RMSEA lower 95% CI 0.0217
## RMSEA upper 95% CI 0.0500
## SRMR 0.0160
## TLI 0.9801
## CFI 0.9920
## AIC 13.9954
## BIC -10.8233
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
##
## $`Younger participants`
##
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj),
## n = sum(tssem1.obj$n), Amatrix = Amatrix, Smatrix = Smatrix,
## Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis,
## intervals.type = intervals.type, mx.algebras = mx.algebras,
## model.name = model.name, suppressWarnings = suppressWarnings,
## silent = silent, run = run)
##
## 95% confidence intervals: z statistic approximation
## Coefficients:
## Estimate Std.Error lbound ubound z value Pr(>|z|)
## AonAlpha 0.747647 0.023880 0.700842 0.794451 31.3081 <2e-16 ***
## ConAlpha 0.911705 0.019864 0.872772 0.950638 45.8968 <2e-16 ***
## EonBeta 0.152561 0.159140 -0.159347 0.464469 0.9587 0.3377
## ESonAlpha 0.677434 0.025864 0.626742 0.728126 26.1924 <2e-16 ***
## IonBeta 3.283875 3.363587 -3.308633 9.876384 0.9763 0.3289
## AlphawithBeta 0.117255 0.125430 -0.128584 0.363094 0.9348 0.3499
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Goodness-of-fit indices:
## Value
## Sample size 838.0000
## Chi-square of target model 145.6167
## DF of target model 4.0000
## p value of target model 0.0000
## Number of constraints imposed on "Smatrix" 0.0000
## DF manually adjusted 0.0000
## Chi-square of independence model 2480.2352
## DF of independence model 10.0000
## RMSEA 0.2057
## RMSEA lower 95% CI 0.1778
## RMSEA upper 95% CI 0.2350
## SRMR 0.1051
## TLI 0.8567
## CFI 0.9427
## AIC 137.6167
## BIC 118.6926
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Setup two plots
layout(t(1:2))
invisible(lapply(fixed2.cluster, plot))
## Setup two plots
layout(t(1:2))
## Plot the first group
plot(fixed2.cluster[[1]])
title("Younger participants")
## Plot the second group
plot(fixed2.cluster[[2]])
title("Older participants")
We may specify RE.type="Diag" with only 10 variances in the diagonal variance component matrix.
The random-effects model is more appropriate.
## method="REM": Random-effects model
random1 <- tssem1(Digman97.data, Digman97.n, method="REM", RE.type="Diag")
summary(random1)
##
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues,
## "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound,
## I2 = I2, model.name = model.name, suppressWarnings = TRUE,
## silent = silent, run = run)
##
## 95% confidence intervals: z statistic approximation
## Coefficients:
## Estimate Std.Error lbound ubound z value Pr(>|z|)
## Intercept1 0.38971902 0.05429255 0.28330757 0.49613047 7.1781 7.068e-13 ***
## Intercept2 0.43265876 0.04145135 0.35141562 0.51390190 10.4377 < 2.2e-16 ***
## Intercept3 0.04945623 0.06071078 -0.06953471 0.16844718 0.8146 0.41529
## Intercept4 0.09603703 0.04478710 0.00825593 0.18381813 2.1443 0.03201 *
## Intercept5 0.42724236 0.03911620 0.35057603 0.50390870 10.9224 < 2.2e-16 ***
## Intercept6 0.11929312 0.04106203 0.03881303 0.19977322 2.9052 0.00367 **
## Intercept7 0.19292424 0.04757961 0.09966991 0.28617856 4.0548 5.018e-05 ***
## Intercept8 0.22690158 0.03240893 0.16338125 0.29042191 7.0012 2.538e-12 ***
## Intercept9 0.18105561 0.04258855 0.09758360 0.26452763 4.2513 2.126e-05 ***
## Intercept10 0.43614966 0.03205960 0.37331400 0.49898532 13.6043 < 2.2e-16 ***
## Tau2_1_1 0.03648370 0.01513054 0.00682839 0.06613902 2.4113 0.01590 *
## Tau2_2_2 0.01963096 0.00859788 0.00277942 0.03648250 2.2832 0.02242 *
## Tau2_3_3 0.04571436 0.01952284 0.00745029 0.08397842 2.3416 0.01920 *
## Tau2_4_4 0.02236119 0.00995082 0.00285794 0.04186445 2.2472 0.02463 *
## Tau2_5_5 0.01729071 0.00796404 0.00168149 0.03289994 2.1711 0.02992 *
## Tau2_6_6 0.01815481 0.00895896 0.00059557 0.03571405 2.0264 0.04272 *
## Tau2_7_7 0.02604881 0.01130265 0.00389601 0.04820160 2.3047 0.02119 *
## Tau2_8_8 0.00988761 0.00513713 -0.00018098 0.01995619 1.9247 0.05426 .
## Tau2_9_9 0.01988243 0.00895052 0.00233973 0.03742513 2.2214 0.02633 *
## Tau2_10_10 0.01049221 0.00501578 0.00066147 0.02032295 2.0918 0.03645 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Q statistic on the homogeneity of effect sizes: 1220.332
## Degrees of freedom of the Q statistic: 130
## P value of the Q statistic: 0
##
## Heterogeneity indices (based on the estimated Tau2):
## Estimate
## Intercept1: I2 (Q statistic) 0.9326
## Intercept2: I2 (Q statistic) 0.8872
## Intercept3: I2 (Q statistic) 0.9325
## Intercept4: I2 (Q statistic) 0.8703
## Intercept5: I2 (Q statistic) 0.8797
## Intercept6: I2 (Q statistic) 0.8480
## Intercept7: I2 (Q statistic) 0.8887
## Intercept8: I2 (Q statistic) 0.7669
## Intercept9: I2 (Q statistic) 0.8590
## Intercept10: I2 (Q statistic) 0.8193
##
## Number of studies (or clusters): 14
## Number of observed statistics: 140
## Number of estimated parameters: 20
## Degrees of freedom: 120
## -2 log likelihood: -112.4196
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))
## Intercept1 Intercept2 Intercept3 Intercept4 Intercept5 Intercept6 Intercept7 Intercept8 Intercept9 Intercept10
## 0.38971902 0.43265876 0.04945623 0.09603703 0.42724236 0.11929312 0.19292424 0.22690158 0.18105561 0.43614966
## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
vec2symMat(est_fixed, diag=FALSE)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.00000000 0.3897190 0.4326588 0.04945623 0.09603703
## [2,] 0.38971902 1.0000000 0.4272424 0.11929312 0.19292424
## [3,] 0.43265876 0.4272424 1.0000000 0.22690158 0.18105561
## [4,] 0.04945623 0.1192931 0.2269016 1.00000000 0.43614966
## [5,] 0.09603703 0.1929242 0.1810556 0.43614966 1.00000000
The stage 2 analysis is similar to that for the fixed-effects model.
The estimated correlation between the two factors (labeled cor in the printout) is \(\hat{S}[7,6]=0.3777\).
random2 <- tssem2(random1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)
summary(random2)
##
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, Amatrix = Amatrix,
## Smatrix = Smatrix, Fmatrix = Fmatrix, diag.constraints = diag.constraints,
## cor.analysis = cor.analysis, intervals.type = intervals.type,
## mx.algebras = mx.algebras, model.name = model.name, suppressWarnings = suppressWarnings,
## silent = silent, run = run)
##
## 95% confidence intervals: z statistic approximation
## Coefficients:
## Estimate Std.Error lbound ubound z value Pr(>|z|)
## AonAlpha 0.569435 0.052425 0.466684 0.672187 10.8619 < 2.2e-16 ***
## ConAlpha 0.590630 0.052649 0.487439 0.693821 11.2182 < 2.2e-16 ***
## EonBeta 0.679955 0.075723 0.531541 0.828370 8.9795 < 2.2e-16 ***
## ESonAlpha 0.760455 0.061963 0.639009 0.881901 12.2726 < 2.2e-16 ***
## IonBeta 0.641842 0.072459 0.499825 0.783859 8.8580 < 2.2e-16 ***
## AlphawithBeta 0.377691 0.047402 0.284785 0.470596 7.9679 1.554e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Goodness-of-fit indices:
## Value
## Sample size 4496.0000
## Chi-square of target model 7.8204
## DF of target model 4.0000
## p value of target model 0.0984
## Number of constraints imposed on "Smatrix" 0.0000
## DF manually adjusted 0.0000
## Chi-square of independence model 501.6767
## DF of independence model 10.0000
## RMSEA 0.0146
## RMSEA lower 95% CI 0.0000
## RMSEA upper 95% CI 0.0297
## SRMR 0.0436
## TLI 0.9806
## CFI 0.9922
## AIC -0.1796
## BIC -25.8234
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Plot the parameter estimates
plot(random2, color="green")
metaSEM package as Becker94.## Read the first sheet in the Excel file
my.xlsx <- read_excel("Becker94.xlsx")
my.xlsx
## # A tibble: 10 x 6
## Study r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal n gender
## <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Becker (1978) Females 0.47 -0.21 -0.15 74 Females
## 2 Becker (1978) Males 0.28 0.19 0.18 153 Males
## 3 Berry (1957) Females 0.48 0.41 0.26 48 Females
## 4 Berry (1957) Males 0.37 0.4 0.27 55 Males
## 5 Rosenberg (1981) Females 0.42 0.48 0.23 51 Females
## 6 Rosenberg (1981) Males 0.41 0.74 0.44 18 Males
## 7 Weiner A (1984) Females 0.26 0.72 0.36 27 Females
## 8 Weiner A (1984) Males 0.32 0.52 0.1 43 Males
## 9 Weiner B (1984) Females 0.580 0.64 0.4 35 Females
## 10 Weiner B (1984) Males 0.34 0.28 -0.03 34 Males
## Split the data by rows
my.R2 <- split(my.xlsx[, 2:4], seq(nrow(my.xlsx)))
head(my.R2)
## $`1`
## # A tibble: 1 x 3
## r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
## <dbl> <dbl> <dbl>
## 1 0.47 -0.21 -0.15
##
## $`2`
## # A tibble: 1 x 3
## r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
## <dbl> <dbl> <dbl>
## 1 0.28 0.19 0.18
##
## $`3`
## # A tibble: 1 x 3
## r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
## <dbl> <dbl> <dbl>
## 1 0.48 0.41 0.26
##
## $`4`
## # A tibble: 1 x 3
## r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
## <dbl> <dbl> <dbl>
## 1 0.37 0.4 0.27
##
## $`5`
## # A tibble: 1 x 3
## r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
## <dbl> <dbl> <dbl>
## 1 0.42 0.48 0.23
##
## $`6`
## # A tibble: 1 x 3
## r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
## <dbl> <dbl> <dbl>
## 1 0.41 0.74 0.44
## Convert the row correlations into correlation matrices
my.R2 <- lapply(my.R2, function(x) vec2symMat(unlist(x), diag=FALSE))
head(my.R2)
## $`1`
## [,1] [,2] [,3]
## [1,] 1.00 0.47 -0.21
## [2,] 0.47 1.00 -0.15
## [3,] -0.21 -0.15 1.00
##
## $`2`
## [,1] [,2] [,3]
## [1,] 1.00 0.28 0.19
## [2,] 0.28 1.00 0.18
## [3,] 0.19 0.18 1.00
##
## $`3`
## [,1] [,2] [,3]
## [1,] 1.00 0.48 0.41
## [2,] 0.48 1.00 0.26
## [3,] 0.41 0.26 1.00
##
## $`4`
## [,1] [,2] [,3]
## [1,] 1.00 0.37 0.40
## [2,] 0.37 1.00 0.27
## [3,] 0.40 0.27 1.00
##
## $`5`
## [,1] [,2] [,3]
## [1,] 1.00 0.42 0.48
## [2,] 0.42 1.00 0.23
## [3,] 0.48 0.23 1.00
##
## $`6`
## [,1] [,2] [,3]
## [1,] 1.00 0.41 0.74
## [2,] 0.41 1.00 0.44
## [3,] 0.74 0.44 1.00
## Add the labels
var.names <- c("SAT_Math", "Spatial", "SAT_Verbal")
my.R2 <- lapply( my.R2, function(x) { dimnames(x) <- list(var.names, var.names); x} )
## Add the study name
names(my.R2) <- my.xlsx$Study
head(my.R2)
## $`Becker (1978) Females`
## SAT_Math Spatial SAT_Verbal
## SAT_Math 1.00 0.47 -0.21
## Spatial 0.47 1.00 -0.15
## SAT_Verbal -0.21 -0.15 1.00
##
## $`Becker (1978) Males`
## SAT_Math Spatial SAT_Verbal
## SAT_Math 1.00 0.28 0.19
## Spatial 0.28 1.00 0.18
## SAT_Verbal 0.19 0.18 1.00
##
## $`Berry (1957) Females`
## SAT_Math Spatial SAT_Verbal
## SAT_Math 1.00 0.48 0.41
## Spatial 0.48 1.00 0.26
## SAT_Verbal 0.41 0.26 1.00
##
## $`Berry (1957) Males`
## SAT_Math Spatial SAT_Verbal
## SAT_Math 1.00 0.37 0.40
## Spatial 0.37 1.00 0.27
## SAT_Verbal 0.40 0.27 1.00
##
## $`Rosenberg (1981) Females`
## SAT_Math Spatial SAT_Verbal
## SAT_Math 1.00 0.42 0.48
## Spatial 0.42 1.00 0.23
## SAT_Verbal 0.48 0.23 1.00
##
## $`Rosenberg (1981) Males`
## SAT_Math Spatial SAT_Verbal
## SAT_Math 1.00 0.41 0.74
## Spatial 0.41 1.00 0.44
## SAT_Verbal 0.74 0.44 1.00
Math, Verbal, and Spatial.## Regression model
model <- "SAT_Math ~ Spatial + SAT_Verbal
## Correlation between Spatial and Verbal
Spatial ~~ SAT_Verbal
## Fix the variances of the independent variables at 1.0
Spatial ~~ 1*Spatial
SAT_Verbal ~~ 1*SAT_Verbal"
## Plot the model
plot(model, color="yellow")
## Convert the equation into RAM
RAM <- lavaan2RAM(model, obs.variables=c("SAT_Math", "Spatial", "SAT_Verbal"))
RAM
## $A
## SAT_Math Spatial SAT_Verbal
## SAT_Math "0" "0*SAT_MathONSpatial" "0*SAT_MathONSAT_Verbal"
## Spatial "0" "0" "0"
## SAT_Verbal "0" "0" "0"
##
## $S
## SAT_Math Spatial SAT_Verbal
## SAT_Math "0*SAT_MathWITHSAT_Math" "0" "0"
## Spatial "0" "1" "0*SpatialWITHSAT_Verbal"
## SAT_Verbal "0" "0*SpatialWITHSAT_Verbal" "1"
##
## $F
## SAT_Math Spatial SAT_Verbal
## SAT_Math 1 0 0
## Spatial 0 1 0
## SAT_Verbal 0 0 1
##
## $M
## SAT_Math Spatial SAT_Verbal
## 1 0 0 0
## method="FEM": fixed-effects TSSEM
fixed1 <- tssem1(Becker94$data, Becker94$n, method="FEM")
## summary of the findings
summary(fixed1)
##
## Call:
## tssem1FEM(Cov = Cov, n = n, cor.analysis = cor.analysis, model.name = model.name,
## cluster = cluster, suppressWarnings = suppressWarnings, silent = silent,
## run = run)
##
## Coefficients:
## Estimate Std.Error z value Pr(>|z|)
## S[1,2] 0.379961 0.037123 10.2351 < 2.2e-16 ***
## S[1,3] 0.334562 0.039947 8.3751 < 2.2e-16 ***
## S[2,3] 0.176461 0.042334 4.1683 3.069e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Goodness-of-fit indices:
## Value
## Sample size 538.0000
## Chi-square of target model 63.6553
## DF of target model 27.0000
## p value of target model 0.0001
## Chi-square of independence model 207.7894
## DF of independence model 30.0000
## RMSEA 0.1590
## RMSEA lower 95% CI 0.1096
## RMSEA upper 95% CI 0.2117
## SRMR 0.1586
## TLI 0.7709
## CFI 0.7938
## AIC 9.6553
## BIC -106.1169
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
We may specify RE.type="Diag" since there are only ten studies.
SAT (Math) and Spatial ability) and 0.7224 (cor between SAT (Math) and SAT (Verbal)), respectively.The random-effects model is more appropriate.
## method="REM": Random-effects model
random1 <- tssem1(Becker94$data, Becker94$n, method="REM", RE.type="Diag")
summary(random1)
##
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues,
## "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound,
## I2 = I2, model.name = model.name, suppressWarnings = TRUE,
## silent = silent, run = run)
##
## 95% confidence intervals: z statistic approximation
## Coefficients:
## Estimate Std.Error lbound ubound z value Pr(>|z|)
## Intercept1 0.3777491 0.0395030 0.3003246 0.4551736 9.5625 < 2.2e-16 ***
## Intercept2 0.3807843 0.0784956 0.2269357 0.5346329 4.8510 1.228e-06 ***
## Intercept3 0.1704927 0.0513545 0.0698398 0.2711457 3.3199 0.0009004 ***
## Tau2_1_1 0.0005038 0.0042009 -0.0077298 0.0087374 0.1199 0.9045407
## Tau2_2_2 0.0416264 0.0257388 -0.0088206 0.0920734 1.6173 0.1058209
## Tau2_3_3 0.0067540 0.0102792 -0.0133928 0.0269008 0.6571 0.5111466
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Q statistic on the homogeneity of effect sizes: 61.02636
## Degrees of freedom of the Q statistic: 27
## P value of the Q statistic: 0.0001932118
##
## Heterogeneity indices (based on the estimated Tau2):
## Estimate
## Intercept1: I2 (Q statistic) 0.0337
## Intercept2: I2 (Q statistic) 0.7224
## Intercept3: I2 (Q statistic) 0.2676
##
## Number of studies (or clusters): 10
## Number of observed statistics: 30
## Number of estimated parameters: 6
## Degrees of freedom: 24
## -2 log likelihood: -22.61046
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))
## Intercept1 Intercept2 Intercept3
## 0.3777491 0.3807843 0.1704927
## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
vec2symMat(est_fixed, diag=FALSE)
## [,1] [,2] [,3]
## [1,] 1.0000000 0.3777491 0.3807843
## [2,] 0.3777491 1.0000000 0.1704927
## [3,] 0.3807843 0.1704927 1.0000000
random2 <- tssem2(random1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)
summary(random2)
##
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, Amatrix = Amatrix,
## Smatrix = Smatrix, Fmatrix = Fmatrix, diag.constraints = diag.constraints,
## cor.analysis = cor.analysis, intervals.type = intervals.type,
## mx.algebras = mx.algebras, model.name = model.name, suppressWarnings = suppressWarnings,
## silent = silent, run = run)
##
## 95% confidence intervals: z statistic approximation
## Coefficients:
## Estimate Std.Error lbound ubound z value Pr(>|z|)
## SAT_MathONSAT_Verbal 0.325853 0.080000 0.169056 0.482649 4.0732 4.638e-05 ***
## SAT_MathONSpatial 0.322194 0.042415 0.239063 0.405325 7.5963 3.042e-14 ***
## SpatialWITHSAT_Verbal 0.170493 0.051355 0.069840 0.271146 3.3199 0.0009004 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Goodness-of-fit indices:
## Value
## Sample size 538.00
## Chi-square of target model 0.00
## DF of target model 0.00
## p value of target model 0.00
## Number of constraints imposed on "Smatrix" 0.00
## DF manually adjusted 0.00
## Chi-square of independence model 110.81
## DF of independence model 3.00
## RMSEA 0.00
## RMSEA lower 95% CI 0.00
## RMSEA upper 95% CI 0.00
## SRMR 0.00
## TLI -Inf
## CFI 1.00
## AIC 0.00
## BIC 0.00
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## S matrix
mxEval(Smatrix, random2$mx.fit)
## [,1] [,2] [,3]
## [1,] 0.7542121 0.0000000 0.0000000
## [2,] 0.0000000 1.0000000 0.1704927
## [3,] 0.0000000 0.1704927 1.0000000
## R2
mxEval(1-Smatrix, random2$mx.fit)[1,1]
## [1] 0.2457879
## Plot the model
plot(random2, color="green")
metaSEM package as Hunter83.## Sheets including correlation matrices
sheet <- as.character(1:14)
my.R2 <- list()
for (i in sheet) {
my.R2[[i]] <- readxl::read_excel("Hunter83.xlsx", sheet=i)
}
## Read study names and sample sizes in sheet "0"
my.study <- readxl::read_excel("Hunter83.xlsx", sheet="0")
my.R2 <- lapply(my.R2, function(x) {x <- unlist(x)
x <- matrix(x, ncol=4)
x <- vechs(x)
vec2symMat(x, diag=FALSE)})
var.names <- c("Ability", "Knowledge", "Work sample", "Supervisor")
my.R2 <- lapply(my.R2, function(x) { dimnames(x) <- list(var.names, var.names); x})
names(my.R2) <- my.study$Study
my.R2[1:2]
## $`Campbell et al. (1973)`
## Ability Knowledge Work sample Supervisor
## Ability 1.00 0.65 0.48 0.33
## Knowledge 0.65 1.00 0.52 0.40
## Work sample 0.48 0.52 1.00 0.23
## Supervisor 0.33 0.40 0.23 1.00
##
## $`Corts et al. (1977)`
## Ability Knowledge Work sample Supervisor
## Ability 1.00 0.53 0.50 0.03
## Knowledge 0.53 1.00 0.49 0.04
## Work sample 0.50 0.49 1.00 0.18
## Supervisor 0.03 0.04 0.18 1.00
model <- "## Regression paths
Job_knowledge ~ A2J*Ability
Work_sample ~ A2W*Ability + J2W*Job_knowledge
Supervisor ~ J2S*Job_knowledge + W2S*Work_sample
## Fix the variance of Ability at 1
Ability ~~ 1*Ability
## Label the error variances of the dependent variables
Job_knowledge ~~ Var_e_J*Job_knowledge
Work_sample ~~ Var_e_W*Work_sample
Supervisor ~~ Var_e_S*Supervisor"
plot(model, color="yellow", layout="spring")
RAM <- lavaan2RAM(model, obs.variables=c("Ability","Job_knowledge",
"Work_sample","Supervisor"))
RAM
## $A
## Ability Job_knowledge Work_sample Supervisor
## Ability "0" "0" "0" "0"
## Job_knowledge "0*A2J" "0" "0" "0"
## Work_sample "0*A2W" "0*J2W" "0" "0"
## Supervisor "0" "0*J2S" "0*W2S" "0"
##
## $S
## Ability Job_knowledge Work_sample Supervisor
## Ability "1" "0" "0" "0"
## Job_knowledge "0" "0*Var_e_J" "0" "0"
## Work_sample "0" "0" "0*Var_e_W" "0"
## Supervisor "0" "0" "0" "0*Var_e_S"
##
## $F
## Ability Job_knowledge Work_sample Supervisor
## Ability 1 0 0 0
## Job_knowledge 0 1 0 0
## Work_sample 0 0 1 0
## Supervisor 0 0 0 1
##
## $M
## Ability Job_knowledge Work_sample Supervisor
## 1 0 0 0 0
We may specify RE.type="Diag" since there are only 14 studies.
The random-effects model is more appropriate.
## method="REM": Random-effects model
random1 <- tssem1(Hunter83$data, Hunter83$n, method="REM", RE.type="Diag")
summary(random1)
##
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues,
## "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound,
## I2 = I2, model.name = model.name, suppressWarnings = TRUE,
## silent = silent, run = run)
##
## 95% confidence intervals: z statistic approximation
## Coefficients:
## Estimate Std.Error lbound ubound z value Pr(>|z|)
## Intercept1 5.0045e-01 2.8249e-02 4.4508e-01 5.5581e-01 17.7157 < 2.2e-16 ***
## Intercept2 4.3124e-01 2.1560e-02 3.8898e-01 4.7350e-01 20.0020 < 2.2e-16 ***
## Intercept3 2.0615e-01 2.5700e-02 1.5578e-01 2.5652e-01 8.0214 1.110e-15 ***
## Intercept4 5.2074e-01 3.1799e-02 4.5842e-01 5.8307e-01 16.3758 < 2.2e-16 ***
## Intercept5 2.5781e-01 3.5914e-02 1.8742e-01 3.2820e-01 7.1785 7.045e-13 ***
## Intercept6 2.3471e-01 1.7505e-02 2.0040e-01 2.6902e-01 13.4085 < 2.2e-16 ***
## Tau2_1_1 6.2569e-03 3.6590e-03 -9.1449e-04 1.3428e-02 1.7100 0.08726 .
## Tau2_2_2 2.2234e-03 2.1500e-03 -1.9905e-03 6.4373e-03 1.0342 0.30106
## Tau2_3_3 4.3637e-03 3.0130e-03 -1.5417e-03 1.0269e-02 1.4483 0.14754
## Tau2_4_4 7.5512e-03 4.6477e-03 -1.5581e-03 1.6661e-02 1.6247 0.10422
## Tau2_5_5 1.0692e-02 6.0658e-03 -1.1966e-03 2.2581e-02 1.7627 0.07795 .
## Tau2_6_6 1.0005e-10 1.7130e-03 -3.3575e-03 3.3575e-03 0.0000 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Q statistic on the homogeneity of effect sizes: 205.225
## Degrees of freedom of the Q statistic: 60
## P value of the Q statistic: 0
##
## Heterogeneity indices (based on the estimated Tau2):
## Estimate
## Intercept1: I2 (Q statistic) 0.7209
## Intercept2: I2 (Q statistic) 0.4475
## Intercept3: I2 (Q statistic) 0.5671
## Intercept4: I2 (Q statistic) 0.7459
## Intercept5: I2 (Q statistic) 0.7691
## Intercept6: I2 (Q statistic) 0.0000
##
## Number of studies (or clusters): 14
## Number of observed statistics: 66
## Number of estimated parameters: 12
## Degrees of freedom: 54
## -2 log likelihood: -130.2043
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Occasionally, we may encounter some error messages.
## We may try to remove the error messages by rerunning the model.
random1 <- rerun(random1)
##
## Beginning initial fit attempt
##
## Solution found
##
## Solution found! Final fit=-130.20428 (started at -130.20428) (1 attempt(s): 1 valid, 0 errors)
summary(random1)
##
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues,
## "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound,
## I2 = I2, model.name = model.name, suppressWarnings = TRUE,
## silent = silent, run = run)
##
## 95% confidence intervals: z statistic approximation
## Coefficients:
## Estimate Std.Error lbound ubound z value Pr(>|z|)
## Intercept1 5.0045e-01 2.8249e-02 4.4508e-01 5.5581e-01 17.7157 < 2.2e-16 ***
## Intercept2 4.3124e-01 2.1560e-02 3.8898e-01 4.7350e-01 20.0020 < 2.2e-16 ***
## Intercept3 2.0615e-01 2.5700e-02 1.5578e-01 2.5652e-01 8.0214 1.110e-15 ***
## Intercept4 5.2074e-01 3.1799e-02 4.5842e-01 5.8307e-01 16.3758 < 2.2e-16 ***
## Intercept5 2.5781e-01 3.5914e-02 1.8742e-01 3.2820e-01 7.1785 7.045e-13 ***
## Intercept6 2.3471e-01 1.7505e-02 2.0040e-01 2.6902e-01 13.4085 < 2.2e-16 ***
## Tau2_1_1 6.2569e-03 3.6590e-03 -9.1449e-04 1.3428e-02 1.7100 0.08726 .
## Tau2_2_2 2.2234e-03 2.1500e-03 -1.9905e-03 6.4373e-03 1.0342 0.30106
## Tau2_3_3 4.3637e-03 3.0130e-03 -1.5417e-03 1.0269e-02 1.4483 0.14754
## Tau2_4_4 7.5512e-03 4.6477e-03 -1.5581e-03 1.6661e-02 1.6247 0.10422
## Tau2_5_5 1.0692e-02 6.0658e-03 -1.1966e-03 2.2581e-02 1.7627 0.07795 .
## Tau2_6_6 1.0005e-10 1.7130e-03 -3.3575e-03 3.3575e-03 0.0000 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Q statistic on the homogeneity of effect sizes: 205.225
## Degrees of freedom of the Q statistic: 60
## P value of the Q statistic: 0
##
## Heterogeneity indices (based on the estimated Tau2):
## Estimate
## Intercept1: I2 (Q statistic) 0.7209
## Intercept2: I2 (Q statistic) 0.4475
## Intercept3: I2 (Q statistic) 0.5671
## Intercept4: I2 (Q statistic) 0.7459
## Intercept5: I2 (Q statistic) 0.7691
## Intercept6: I2 (Q statistic) 0.0000
##
## Number of studies (or clusters): 14
## Number of observed statistics: 66
## Number of estimated parameters: 12
## Degrees of freedom: 54
## -2 log likelihood: -130.2043
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))
## Intercept1 Intercept2 Intercept3 Intercept4 Intercept5 Intercept6
## 0.5004463 0.4312399 0.2061453 0.5207418 0.2578125 0.2347090
## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
vec2symMat(est_fixed, diag=FALSE)
## [,1] [,2] [,3] [,4]
## [1,] 1.0000000 0.5004463 0.4312399 0.2061453
## [2,] 0.5004463 1.0000000 0.5207418 0.2578125
## [3,] 0.4312399 0.5207418 1.0000000 0.2347090
## [4,] 0.2061453 0.2578125 0.2347090 1.0000000
Ind = A2J*J2S+A2J*J2W*W2S+A2W*W2Sintervals.type="LB".random2 <- tssem2(random1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F,
intervals.type="LB",
mx.algebras=
list( ind=mxAlgebra(A2J*J2S+A2J*J2W*W2S+A2W*W2S, name="ind") ) )
summary(random2)
##
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, Amatrix = Amatrix,
## Smatrix = Smatrix, Fmatrix = Fmatrix, diag.constraints = diag.constraints,
## cor.analysis = cor.analysis, intervals.type = intervals.type,
## mx.algebras = mx.algebras, model.name = model.name, suppressWarnings = suppressWarnings,
## silent = silent, run = run)
##
## 95% confidence intervals: Likelihood-based statistic
## Coefficients:
## Estimate Std.Error lbound ubound z value Pr(>|z|)
## A2J 0.510783 NA 0.456554 0.564987 NA NA
## J2S 0.224033 NA 0.136802 0.311682 NA NA
## W2S 0.121026 NA 0.056817 0.181657 NA NA
## A2W 0.230913 NA 0.157756 0.299775 NA NA
## J2W 0.397800 NA 0.311509 0.484687 NA NA
##
## mxAlgebras objects (and their 95% likelihood-based CIs):
## lbound Estimate ubound
## ind[1,1] 0.1374292 0.1669702 0.1983626
##
## Goodness-of-fit indices:
## Value
## Sample size 3975.0000
## Chi-square of target model 3.5748
## DF of target model 1.0000
## p value of target model 0.0587
## Number of constraints imposed on "Smatrix" 0.0000
## DF manually adjusted 0.0000
## Chi-square of independence model 975.2843
## DF of independence model 6.0000
## RMSEA 0.0255
## RMSEA lower 95% CI 0.0000
## RMSEA upper 95% CI 0.0561
## SRMR 0.0204
## TLI 0.9841
## CFI 0.9973
## AIC 1.5748
## BIC -4.7129
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Display the A matrix
mxEval(Amatrix, random2$mx.fit)
## Ability Job_knowledge Work_sample Supervisor
## Ability 0.0000000 0.0000000 0.0000000 0
## Job_knowledge 0.5107833 0.0000000 0.0000000 0
## Work_sample 0.2309135 0.3978003 0.0000000 0
## Supervisor 0.0000000 0.2240334 0.1210258 0
## Display the S matrix
mxEval(Smatrix, random2$mx.fit)
## [,1] [,2] [,3] [,4]
## [1,] 1 0.0000000 0.0000000 0.000000
## [2,] 0 0.7391004 0.0000000 0.000000
## [3,] 0 0.0000000 0.6945954 0.000000
## [4,] 0 0.0000000 0.0000000 0.907194
plot(random2, color="yellow", layout="spring")
National Research Council (1992). Combining information: Statistical issues and opportunities for research. Washington, D.C.: National Academy Press.↩
Brown, S. P., & Stayman, D. M. (1992). Antecedents and consequences of attitude toward the ad: A meta-analysis. Journal of Consumer Research, 19, 34-51.↩
Premack, S. L., & Hunter, J. E. (1988). Individual unionization decisions. Psychological Bulletin, 103, 223-234.↩
Norton, S., Cosco, T., Doyle, F., Done, J., & Sacker, A. (2013). The Hospital Anxiety and Depression Scale: A meta confirmatory factor analysis. Journal of Psychosomatic Research, 74(1), 74-81. http://doi.org/10.1016/j.jpsychores.2012.10.010↩
Murayama, K., & Elliot, A. J. (2012). The competition-performance relation: A meta-analytic review and test of the opposing processes model of competition and performance. Psychological Bulletin, 138(6), 1035-1070. http://doi.org/10.1037/a0028324↩
Viswesvaran, C., & Ones, D. S. (1995). Theory testing: Combining psychometric meta-analysis and structural equations modeling. Personnel Psychology, 48(4), 865-885.↩
Becker, B. J. (1992). Using results from replicated studies to estimate linear models. Journal of Educational Statistics, 17(4), 341-362.↩
Cheung, M. W.-L. (2014). Fixed- and random-effects meta-analytic structural equation modeling: Examples and analyses in R. Behavior Research Methods, 46(1), 29-40. http://doi.org/10.3758/s13428-013-0361-y↩
Cheung, M. W.-L., & Chan, W. (2005). Meta-analytic structural equation modeling: A two-stage approach. Psychological Methods, 10(1), 40-64. http://doi.org/10.1037/1082-989X.10.1.40↩
Cheung, M. W.-L. (2018). Some reflections on combining meta-analysis and structural equation modeling. Research Synthesis Methods, 0(0). https://doi.org/10.1002/jrsm.1321↩
Cudeck, R. (1989). Analysis of correlation matrices using covariance structure models. Psychological Bulletin, 105(2), 317-327.↩
Joreskog, K. G., & Sorbom, D. (1996). LISREL 8: A user-s reference guide. Chicago, IL: Scientific Software International, Inc.↩
Muthen, B., Kaplan, D., & Hollis, M. (1987). On structural equation modeling with data that are not missing completely at random. Psychometrika, 52(3), 431-462. http://doi.org/10.1007/BF02294365↩
Cheung, M. W.-L., & Chan, W. (2009). A two-stage approach to synthesizing covariance matrices in meta-analytic structural equation modeling. Structural Equation Modeling: A Multidisciplinary Journal, 16(1), 28-53. http://doi.org/10.1080/10705510802561295↩
Bentler, P. M., & Savalei, V. (2010). Analysis of correlation structures: Current status and open problems. In S. Kolenikov, D. Steinley, & L. Thombs (Eds.), Statistics in the Social Sciences (pp. 1-36). New Jersey: John Wiley & Sons, Inc.↩
Cheung, M. W.-L., & Chan, W. (2005). Classifying correlation matrices into relatively homogeneous subgroups: A cluster analytic approach. Educational and Psychological Measurement, 65(6), 954–979. https://doi.org/10.1177/0013164404273946↩
Jak, S., & Cheung, M. W.-L. (2018). Testing moderator hypotheses in meta-analytic structural equation modeling using subgroup analysis. Behavior Research Methods, 50(4), 1359–1373. https://doi.org/10.3758/s13428-018-1046-3↩
Rosenthal, R. (1991). Meta-analytic procedures for social research (Revised). Newbury Park, [Calif.]: Sage.↩
Schmidt, F. L., & Hunter, J. E. (2015). Methods of meta-analysis: Correcting error and bias in research findings (3rd ed.). Thousand Oaks, CA: Sage.↩
Michel, J. S., Viswesvaran, C., & Thomas, J. (2011). Conclusions from meta-analytic structural equation models generally do not change due to corrections for study artifacts. Research Synthesis Methods, 2(3), 174-187. http://doi.org/10.1002/jrsm.47↩
McArdle, J. J., & McDonald, R. P. (1984). Some algebraic properties of the Reticular Action Model for moment structures. British Journal of Mathematical and Statistical Psychology, 37(2), 234-251. http://doi.org/10.1111/j.2044-8317.1984.tb00802.x↩
Browne, M. W., & Cudeck, R. (1993). Alternative ways of assessing model fit. In K. A. Bollen & J. S. Long (Eds.), Testing Structural Equation Models (pp. 136-162). Newbury Park, CA: Sage.↩
Digman, J. M. (1997). Higher-order factors of the Big Five. Journal of Personality and Social Psychology, 73(6), 1246-1256. http://doi.org/10.1037/0022-3514.73.6.1246↩
Epskamp, S (2017). semPlot: Path diagrams and visual analysis of various SEM packages’ output. R package version 1.1.0. http://CRAN.R-project.org/package=semPlot↩
Becker, B. J., & Schram, C. M. (1994). Examining explanatory models through research synthesis. In H. Cooper & L. V. Hedges (Eds.), The handbook of research synthesis (pp. 357-381). New York: Russell Sage Foundation.↩